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,^,! Abstract 

The multilayer multiconfiguration time-dependent Hartree method is employed to study vibra- 

■ tionally coupled charge transport in models of single molecule junctions. To increase the efficiency 

cn ■ 

■ of the simulation method, a representation of the Hamiltonian in terms of the scattering states 

cn ■ 

(N 



of the underlying electronic Hamiltonian is used. It is found that with an appropriate choice of 



I the scattering states the artificial electron correlation present in the original representation of the 

I model is greatly reduced. This allows efficient simulation of the steady-state currents in a wide 

physical parameter space, which is demonstrated by several numerical examples. 
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I. INTRODUCTION 



Charge transport in single-molecule junctions has raised a great deal of interest 
recently.-"— Experimental studies on transport properties of molecular junctions revealed 
a variety of interesting phenomena, e.g., Coulomb blockade,— Kondo effect,™ negative dif- 
ferential resistance,—""— switching and hysteresis.—"— This has stimulated many theoretical 
developments for understanding quantum transport at the molecular scale. Examples of 
approximate methods employed in this regard include the scattering theory,—"— nonequilib- 
rium Green's function (NEGF) approaches,—"— and master equation methods.— i^^"— To the 
other end, a variety of (in principle) numerically exact methods have been developed to ob- 
tain more reliable results for nonequilibrium transport in model systems. These include the 
numerical path integral approach,—"— real-time quantum Monte Carlo simulations,—!^ the 
numerical renormalization group approach,Si, the time-dependent density matrix renormal- 
ization group approach. and the hierarchical equations of motion method.— Our work 
in this direction is the development of the multilayer multiconfiguration time-dependent 
Hartree (ML-MCTDH) theory in the second quantization representation (SQR),— a system- 
atic, numerically exact methodology to study quantum dynamics and quantum transport 
including many-body effects. For a generic model of vibrationally coupled electron trans- 
port, we have demonstrated^^ the importance of treating the vibronic coupling accurately. 
Comparison with approximate methods such as NEGF revealed the necessity of employ- 
ing accurate methods such as the ML-MCTDH-SQR approach, in particular in the strong 
coupling regime. 

The ML-MCTDH-SQR method^ employs a recursive, layered representation of the over- 
all wave function and captures correlation effects through a converged multiconfigurational 
expansion in all the corresponding layers. Like any other methods, the efficiency of the 
ML-MCTDH-SQR theory depends on how the problem is approached, e.g., the choice of 
the an appropriate coordinate system. Consequently, the term "correlation" also depends 
on such a choice. For example, the quantum dynamics of a quadratic Hamiltonian, where 
a full Hessian is present, is highly correlated. On the other hand, when this Hamiltonian 
is transformed to the normal mode representation, there is no correlation at all between 
different degrees of freedom. In this sense we can term the former correlation as "artificial 
correlation" due to the inappropriate representation of the Hamiltonian. Finding an op- 



timal representation is thus crucial for solving the problem efficiently. Unfortunately, this 
is typically not an automatic procedure and requires a clear physical understanding of the 
problem. In most simple analytic reduction of the problem as the example above 

is not available. Instead, one has to transform the Hamiltonian according to the physical 
intuition and known experiences. One of the most frequently used criterion is to ensure that 
this transformation reduces the problem in a certain physical limit. 

In a previous ML-MCTDH-SQR study of vibrationally coupled electron transport,— we 
had found that in some regimes the configuration space required to achieve convergence is 
quite large. We noticed, however, that these regimes are not necessarily characterized by a 
strong correlation in the usual sense. In fact, it was often relatively easy to converge the 
simulation when the electron- nuclear coupling is strong (i.e., easier to capture the "real" 
correlation effect induced by electronic- vibrational coupling), but not otherwise (i.e. more 
difficult to capture artificial electron correlation effects), especially when the source-drain 
bias voltage is high. This has motivated us to examine other representations of the Hamil- 
tonian for studying nonequilibrium charge transport. One promising choice is to employ the 
scattering states of the underlying electronic Hamiltonian. The resulting scattering states 
representation of the model has been used before to describe quantum transport in com- 
bination with different methods.—"— In this paper we discuss our implementation of the 
scattering state representation within the framework of the ML-MCTDH-SQR theory. 

The paper is organized as follows. Section [TTl outlines the original, tight-binding type 
model for describing the vibrationally coupled electron transport. Section [XTTl then discusses 
the scattering state representation of the same problem and our specific choice of the refer- 
ence frame. Section HV] illustrates the performance of our approach by numerical examples 
for a variety of parameter regimes. Section |VT] concludes with a summary. 

II. MODEL FOR VIBRATIONALLY COUPLED ELECTRON TRANSPORT 
A. Model Hamiltonian 

A model that is frequently adopted to study vibrationally coupled electron transport 
through a single-molecule junction is based on tight-binding description for the two metal 
leads. It comprises one discrete electronic state at the molecular bridge, two electronic 



continua describing the left and the right metal leads, respectively, and a distribution of 
harmonic oscillators that models the vibrational modes of the molecular bridge. The Hamil- 
tonian reads 

H = Hel + -ffel-nuc, (2.1a) 

where He\ and -ffci-nuc describe the pure electronic degrees of freedom and the nuclear vi- 
brations with their coupling terms to the electronic part, respectively 



kh kji 
He\—nuc 



\ E(^' + ^"Q") + J2 ^c,Qr (2.1c) 
j j 

Thereby, d'^/d, 4^/ck^, c^^/cfc^ are the fermionic creation/annihilation operators for the 
electronic states on the molecular bridge, the left and the right leads, respectively. The 
corresponding electronic energies i?fc^, E^j^ and the molecule- lead coupling strengths Vdh^, 
Vdkji, are defined through the energy-dependent level width functions 

TL{E) = 27TY\ydkA^6{E-EkJ, Tr{E) = 27r Y\ydkJ^S{E-Ek^). (2.2) 

kh kn 

In principle, the parameters of the model can be obtained for a specific molecular junction 
employing first-principles electronic structure calculations.— In this paper, which focuses on 
the methodology and general transport properties, however, we will use a generic parame- 
terization. Employing a tight-binding model, the function T{E) is given as 



T{E) = <^ V '^^ I I _ lA^ei ^ ^2.3a) 

\E\>2\P,\ 

Tl{E) = T{E - fiL), TniE) = T{E - ^ir), (2.3b) 

where Pe and ae are nearest-neighbor couplings between two lead sites and between the lead 
and the bridge state, respectively. I.e., the width functions for the left and the right leads 
are obtained by shifting T{E) relative to the chemical potentials of the corresponding leads. 
We consider a simple model of two identical leads, in which the chemical potentials are given 
by 

I^L/R = Ef± r/2, (2.4) 



where V is the source-drain bias voltage and Ef the Fermi energy of the leads. Since only 
the difference Ed — Ef is physically relevant, we set -E/ = 0. 

Similarly, the frequencies Uj and electronic-nuclear coupling constants Cj of the vibrational 
modes of the molecular junctions are modeled by a spectral density function^i^ 




In this paper, the spectral density is chosen in Ohmic form with an exponential cutoff 

Joico) = |aa;e--/'^=, (2.6) 

where a is the dimensionless Kondo parameter. 

Both the electronic and the vibrational continua can be discretized by choosing a density 
of states Pe{E) and a density of frequencies p{ijj) such that^"— 

[ " dE p,{E) = k, |FrffcP = /™, , k=l,...,N,, (2.7a) 

Jo ZnpeiEk) 

rdup{u)=3, c|^2Jo(c^^ j = l,...,iV,. (2.7b) 

where is the number of electronic states (for a single spin/single lead) and is the 
number of bath modes in the simulation. In this work, we choose a constant Pe{E), i.e., an 
equidistant discretization of the interval [— 2/3e, 2/3e], to discretize the electronic continuum. 
For the vibrational bath, p{u) is chosen as 

p(uj) = i!^e-'^/-^ (2.8) 

OJc 

Within a given time scale the numbers of electronic states and bath modes are systematically 
increased to reach converged results for the quantum dynamics in the extended condensed 
phase system. 



B. Initial state and calculation of current 



In this paper the observable of interest for studying transport through molecular junctions 
is the current for a given source-drain bias voltage, given by (we use atomic units where 
h = e = l) 

hit) = = -^tr {pe^^H[H, iV^]e-^^*} , (2.9a) 



lR{t) = = ^tr {pe^^H[H, Nn]e-^^'} . (2.9b) 

Here, Ni/Rit) denotes the time-dependent charge in each lead, defined as 

Ndt) = ^tr[pe^^*iVce-^*], (2.10) 

and = '^f.^ (^t^^^k^ is the occupation number operator for the electrons in each lead 
{( = L,R). For Hamiltonian (12. ip the explicit expression for the current operator is given 



as 



ic^i[H,N^]=iY,Vdk,{d+Ck,-ct^d), C = L,R. (2.11) 

In the expressions above, p denotes the initial density matrix representing a grand-canonical 
ensemble for each lead and a certain preparation for the bridge state 

p = /5° exp [-/3(iJo - f^LNL - PrNr)] , (2.12a) 

Ho = J2 Ek,nk„. + Yl + ^nuc- (2.12b) 

kL,a kR,(7 

Here is the initial reduced density matrix for the bridge state, which is usually chosen as 
a pure state representing an occupied or an empty bridge state, and H^^^ defines the initial 
bath equilibrium distribution. 

Various initial states can be considered. For example, one may choose an initially unoc- 
cupied bridge state and the nuclear degrees of freedom equilibrated with this state, i.e. an 
unshifted bath of oscillators with 

j 

On the other hand, one may also start with a fully occupied bridge state and a bath of 
oscillators in equilibrium with the occupied bridge state 

2" 



mi. = -y 

2 ^ 
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F,^ + c.|(g, + 2-| 



(2.14) 



Other initial states may also be prepared. The initial state may affect the transient dynamics 
profoundly. The dependence of the steady-state current on the initial density matrix is a 
more complex issue. Recent investigations for a model without electron-electron interaction 
seem to indicate that different initial states may lead to different (quasi)steady states,—"-^ 



although this has been debated.— In this paper, we hmit our discussion to physical regimes 
where the steady state current does not depend on the initial state. When comparing the 
time-dependent current from the scattering state representation with that from the original 
Hamiltonian discussed above, we typically choose initial conditions that are close to the final 
steady state, e.g., an unoccupied initial bridge state if its energy is higher than the Fermi 
level of the leads and an occupied bridge state otherwise. We also use the average current 

J(t) = i[/^(t) + /i(t)], (2.15) 

for such a purpose, which minimizes transient effects. 

In our simulations the continuous set of electronic states of the leads is represented by a 
finite number of states. The number of states required to properly describe the continuum 
limit depends on the time t. The situation is thus similar to that of a quantum reactive 
scattering calculation in the presence of a scattering continuum, where, with a finite number 
of basis functions, an appropriate absorbing boundary condition is added to mimic the 
correct outgoing Green's function.—"— Employing the same strategy for the present problem, 
the regularized electric current is given by 

r^s= lim r dt^^e-^^K (2.16) 

The regularization parameter rj is similar (though not identical) to the formal convergence 
parameter in the definition of the Green's function in terms of the time evolution operator 

G{E+) = \im {-i) dte'^^+''^-"^K (2.17) 

'7^0+ Jo 

In numerical calculations, rj is chosen in a similar way as the absorbing potential used in 
quantum scattering calculations.—"— In particular, the parameter rj has to be large enough to 
accelerate the convergence but still sufficiently small in order not to affect the correct result. 
While in the reactive scattering calculation rj is often chosen to be coordinate dependent, in 
our simulation rj is chosen to be time dependent 



i^{t) = { ^ \ (2.18) 




Here r^o is a damping constant, r is a cutoff time beyond which a steady state charge flow 
is approximately reached. As the number of electronic states increases, one may choose a 



weaker damping strength rjo and/or longer cutoff time r. The former approaches zero and 
the latter approaches infinity for an infinite number of states. For the systems considered 
in this work, convergence can be reached with a typical r = 30-80 fs (a smaller r for less 
number of states) and I/tjq = 3-10 fs. 



III. SCATTERING STATE REPRESENTATION OF THE MODEL 

In previous work^i^ we have used the above Hamiltonian with the ML-MCTDH-SQR 
theory for studying vibrationally coupled electron transport through molecular junctions. 
It provided important benchmark results to guide further development of the approximate 
theories. However, it was also found that upon increasing the bias voltage the number of 
configurations required in the ML-MCTDH representation of the wave function to achieve 
convergence becomes very large. This is so even without including the contribution from 
the nuclear bath, Eq. (12. Id) , i.e. for the so-called noninteracting model. In fact, in many 
cases, including the vibrational coupling accelerates the convergence of the ML-MCTDH- 
SQR simulation within a given time scale (see below for a physical explanation). 

Such a performance is not surprising if one considers the electronic Hamiltonian in 
Eq. ( I2.1b|) . Similar to the example of a quadratic Hamiltonian with full Hessian, all the 
electronic degrees of freedom in Eq. fl2.1bl) may be strongly coupled, thus creating signifi- 
cant artificial correlations in the quantum dynamics of the model. Including the vibronic 
contribution in Eq. (I2.1cp may suppress such artificial correlations, though it also introduces 
vibrationally induced 'true' electronic correlation effect. Within the current form of the 
Hamiltonian, these artificial correlations are real and cannot be removed in a variational 
calculation. 

To circumvent this problem and to reduce the artificial correlation, we take advantage of 
the fact that the Green's function and also the steady state current for the noninteracting 
model defined by Eq. (]2.1b|) can be obtained analytically. Thus, the electronic Hamiltonian 
Eq. ( I2.1b[) can be represented in the basis of its eigenstates, which are scattering states. 



The resulting scattering states representation of the model has been used before to describe 
quantum transport in combination with different methods.—"— 



A. General formulation 



In the scattering state representation^"— a set of scattering state operators are introduced 
to diagonalize the Hamiltonian fl2.1b|) . These operators 7^^ satisfy the Lippmann-Schwinger 
equation = L, R) 



— Cfc^ + 

or the equivalent Heisenberg equation 



1 



(3.1a) 



Here = Hq + V, Lq = [Ho,] is a Liouville operator, t] is an infinitesimal, and 

Ho = Ed(td + EkcctcCkc, (3.2a) 



(3.2b) 



The scattering state operators 7^^/7a;c fulfill the usual fermionic commutation relations. 
They can be expanded in terms of and c^^, and the formal solution to the Lippmann- 
Schwinger equation is given by 



(3.3) 



where g%Ek(^) is the retarded Green's function for the bridge state 

1 



9d{Ek() 



Ekc + irj - Ed-T.{Ekc)' 



nEkd = E 



Ekc + iV- Ek'c 



(3.4a) 



(3.4b) 



For the tight-binding parameterization, Eq. (12. 3p . the self energy Ti{E) = Hi^E) + 'Lr{E) 
has the following analytic form 



T.dE) = --TdE) + /\dE), C = L,R, 



(3.5) 



where r^(£') is given by Eq. (12. 3p . and 



A^E) 



2/3; 



iE-f^c) 
(E - /.c) T -{E- ^cf 



1^-/1^1 < 2/3, 
±{E- /i^) > 2/3e 



(3.6) 



In the following, we assume that the purely electronic problem defined by the Hamiltonian 
Eq. fl2.1b[) does not possess bound states, so that the set of the scattering states provide 
a complete basis. As a result, the electronic Hamiltonian is diagonal in the scattering 
representation and has the form 

Hel = ^kai^K = ^kLlt^kr, + ^knlta^kn ■ (3-7) 

For the vibronic part of the Hamiltonian in Eq. f l2.1cp and the current operator in Eq. (12.111) . 
one applies the inverse relation of Eq. (13. 3 p 

^^ = E^£K(^'^c)^'^'^c]* (3.8a) 



and 



9d\^k'C' 



(3.8b) 



which, when substituting into Eq. (12. Ic^ and Eq. ( 12. lip , yields the necessary expressions for 
implementation. It is noted that both the number operator d~^d and the current operator 
are more complex in the scattering state representation. This is the price for reducing the 
artificial correlation effect in the quantum dynamics simulation. 

Finally, the initial density matrix, which will be used in the current simulation based on 
the scattering state representation, is given by 

-/3(Hi-y) 

P= :—, (3.9) 

^ tr[e-/3(^i-^)] ^ 
where Hi is the Hamiltonian without the vibronic coupling 

Hi = Y E^Lltlk, + Y E^nltn^k, + ^ic (3.10) 
kh kji 



and Y is the bias operator 



^ = l^LY^t^kL+ f^nY^k^^f'ri- (3-11) 

kh kn 



The density matrix in Eq. (13. 9p can be factorized 



P — PelPnucb 



(3.12) 



into an electronic part 



g-/3(Hei-y) 



and a nuclear part 

-AH" 

Pnucl = , . .^0 1 - (3-14) 

The electronic density matrix pei describes the steady state of the purely electronic model 
(without electron- vibrational coupling), i.e. 

~ . - f p-/3(ifci-V') 

lim tr{poe'^"'*Ae-'^='n = tr <^ ^ A }> , (3.15) 

holds for a given density matrix po and electronic operator A.—^ In particular, the steady 
state current for vanishing electronic-vibrational coupling is given by 

-/3(H,i-y) ^ 

lim = Ttr ^ ^ . P- 16) 

Assuming a unique steady state, both initial states, Eqs. f l3.9p and (12.121) . will result in the 
same steady-state current, however, the transient dynamics will be different. 



B. Optimization of the transport calculation including vibronic coupling 

It is straightforward to implement the above Hamiltonian in the scattering state repre- 
sentation for the ML-MCTDH-SQR simulation of the current. Without vibronic coupling, 
the model corresponds to a non-interacting system. We have verified that with a sufficient 
number of states used in the discretization of the electronic continua of the leads, the steady 
state current agrees with the Landauer formula, which is exact in this noninteracting case. 
Moreover, the ML-MCTDH-SQR only needs one configuration to obtain the numerically 
exact result, as it should be. 

Including vibronic coupling, this is no longer the case because the systems represents 
a true many-body problem. In this case, we have found that the direct application of 
the scattering state representation as introduced above may result in serious numerical 
problems except for a range of positive energies for the original bridge state {E^ > Ef) where 
limited success is achieved. In particular, when E^ approaches zero or becomes negative, 
the simulated steady-state current sometimes becomes negative even with a relatively large 



number of scattering states and configurations in the ML-MCTDH-SQR calculation, which 
is clearly unphysical. 

This numerical problem may be understood from the technical aspects of applying a fi- 
nite basis set representation to study dynamics. In principle, an infinite basis expansion can 
solve the problem exactly, but in practice one wants to keep the number of basis functions 
as small as possible to make the simulation feasible. To achieve this, it is advantageous to 
keep the (complete) basis sets as close as possible to the eigenfunctions of the overall Hamil- 
tonian. If this is not done appropriately, one may need a prohibitively large number of basis 
functions (i.e., an "infinite" basis in a practical sense) to obtain numerical convergence, such 
as what we have observed here. In the discretized scattering state representation, the over- 
all current comes from the contribution of all the scattering channels. The weight of each 
channel/scattering state is determined properly from the solution to the pure electronic 
Hamiltonian. However, when including the vibrational degrees of freedom in the Hamil- 
tonian and electron-vibrational coupling, such a basis set may be far from optimal. The 
channels that contribute the most to the original, non-interacting electronic system may not 
be important at all to the overall true many-body transport problem. On the other hand, 
other channels that were insignificant previously may now become important. 

According to our investigations, this numerical problem can be circumvented by noticing 
the fact that the scattering state representation is not unique. In particular, there is an 
arbitrariness in defining the bridge state energy, e.g., one can rewrite Eq. (12. ip as 

^ei = {Ea - AE)d+d + EkA.Ck, + EkAn^kn (3.17a) 

kh ku 

i^ei-nuc = I 5Z(P/ + io]Q]) + d+d{AE + Y 2c,Q,). (3.17b) 
j j 
With a different choice of AE, the scattering state operators are different and so will be the 

performance of the vibronic transport calculation. Our numerical tests show that choosing 

specifically the value of AE, which is obtained by completing the square in the vibronic 



part, i.e. 



kL kji 



ka^ka 



(3.18a) 



H, 



cl— nuc 



(3.18b) 



corresponding to Ai? = Ylij results in a particularly good performance of the trans- 

port calculations. It is recognized that for this choice of AE, the correction to the bridge 
state energy, X]j 20^/^1, is the reorganization energy (polaron shift) for electron transfer 
between the bridge and the lead states and thus this selection also has a clear physical 
meaning. 

With the choice of the scattering state operators that diagonalize the H^i in Eq. I3.18a[ 
the effective Hamiltonian for the overall vibronic model performs well in all regimes we have 
tested. The numerical examples presented below will illustrate this point. 



IV. SIMULATION OF TRANSPORT IN THE SCATTERING STATE REPRE- 
SENTATION USING THE ML-MCTDH-SQR METHOD 

To simulate transport in the model described above, we employ the Multilayer Mul- 
ticonfiguration Time-Dependent Hartree Theory in Second Quantization Representation 
(ML-MCTDH-SQR),— which allows a numerically exact treatment of the many-body prob- 
lem. The method has been described in detail elsewhere.— It is based on the ML-MCTDH 
theory^, which is a rigorous variational method to propagate wave packets in complex sys- 
tems with many degrees of freedom. In the ML-MCTDH theory^ the wave function is 
represented by a recursive, layered expansion, 

h 32 jp K=i 

M(K,q) 

i<"^w) = EE- E n i^^f'W)' 

ai a2 "M{K,q) 7=1 



where Aj-^^j^^^j^it), B^^f^^ j^^^^^it) , Caia2...aM{K,q)i't) and so on are the expansion coefficients for 
the ffist, second, third, layers, respectively; \ip^^'^J [t)) , \vll^''^\t)) , \C,a^'''{t)), are the 
"single particle" functions (SPFs) for the ffist, second, third, layers. In Eq. fl4.1al) . p 
denotes the number of single particle (SP) groups/subspaces for the ffist layer. Similarly, 
Q{k) in Eq. (]4.1bp is the number of SP groups for the second layer that belongs to the 
Kth SP group in the ffist layer, i.e., there are a total of X]k=i Q(^) second layer SP groups. 
Continuing along the multilayer hierarchy, M{n, q) in Eq. (14. Id) is the number of SP groups 
for the third layer that belongs to the gth SP group of the second layer and the /tth SP 
group of the ffist layer, resulting in a total of X]k=i Xl^ii'' ^(/^j q) third layer SP groups. 

Each summation in the expressions above is over the number of SPFs in that single 
particle group and the total combination defines the size of the configurational space in that 
layer, e.g., if there are 4 SP groups and we use 10 SPFs per group, the configurational space 
is 10^ = 10000. If there is no correlation at all between different degrees of freedom, only 
one SPF is required for each group, resulting in one overall configuration. As the correlation 
becomes stronger, more SPFs are required to achieve convergence, which results in a larger 
configuration space. Thus, the number of SPFs is a good indicator of the correlation strength 
and will be used for analysis purposes below. 

V. RESULTS AND DISCUSSION 

In this section, we apply the methodology discussed above to study vibrationally-coupled 
charge transport. In particular, we discuss for selected examples the performance of the 
simulation in the scattering-state representation as compared to the original representation. 
Without specifying otherwise, we will use the tight-binding parameterization for the elec- 
tronic part of the model, as described in Sec. [TTland used equivalently in the scattering state 
representation. Sec. Illli The tight-binding parameters for the function T{E) are ae = 0.2 
eV, /3e = 1 eV, corresponding to a moderate molecule-lead coupling an a bandwidth of 4 eV. 



A. Transport without vibronic coupling 



If the Hamiltonian only contains the pure electronic part, i.e. without electronic vibra- 
tional coupling, it corresponds to the non-interacting transport model. Employing the 
scattering-state representation, in this case one SPF for each SP group, i.e., one overall 
configuration representing the wave function, gives the numerically exact steady-state cur- 
rent. This was verified by comparing the numerical results obtained with one configuration 
of the wave function to the Landauer formula. In contrast, for the original representation 
of the model described in Sec. [Tll there is significant artificial correlation, especially for a 
large source-drain bias voltage, such that a large configurational space is required to achieve 
convergence. 

B. Vibrationally coupled electron transport for > Ef 

Including electronic-vibrational coupling, the model becomes a many-body problem with 
"true" correlation. We first consider a physical regime where, without vibronic coupling, the 
electronic energy of the bridge state is located above the Fermi level of the leads, Ed — Ef = 
0.5eV. For small bias voltages, this corresponds to the non-resonant transport regime. Fig.Ht 
shows the converged time- dependent currents obtained with the original model Hamiltonian 
in Sec. [Ill and with the scattering state Hamiltonian. Thereby, the characteristic frequency 
of the bath has been chosen as Uc = 500 cm~^ and the overall electronic-nuclear coupling 
strength is determined by the reorganization energy, A = 2aUc = 2000 cm~^. The source- 
drain bias voltage is O.IV. Despite the difference in the transient behavior of I{t), which is 
expected due to the difference in the initial states between the two model representations, the 
stationary values of the current agree within 1% relative difference. This small numerical 
error is remarkable considering the two drastically different representations of the model 
and the wave function. Fig. [Do shows the convergence behavior of I{t) versus the number 
of SPFs for each electronic SP group, where 6 SPFs for each nuclear SP group provides 
a converged solution. It is seen that the result obtained with 8 SPFs per each electronic 
SP group agrees perfectly with a larger 12 SPFs simulation. Even 6 SPFs per electronic 
SP group gives a result that is within 10% relative error of the converged answer. On the 
other hand, to achieve convergence for the original model requires 40 SPFs per electronic 



SP group to generate the result in Fig. [T^. 

This shows that by using the Hamiltonian in the scattering state representation much of 
the artificial correlation has been removed. This results in a significantly smaller configura- 
tion space than that required for the original model Hamiltonian. A comparison of results 
obtained with different number of SPFs also reveals the extent of the "true" correlation 
effect, which is helpful for analyzing the physics of the problem. 

It is also interesting to consider the result obtained with 1 SPF, i.e., only one overall 
configuration. Within the scattering state representation, this corresponds to the Hartree- 
Fock approximation for the present vibronic model. Fig. [T]d shows that in this case: (i) there 
is no transient dynamics in the simulated I{t) and (ii) the stationary current is approximately 
three times as large as the true result. Further inspection shows that the one-configuration 
result agrees with a pure electronic model (Landauer formula), where the bridge state energy 
is shifted by the bath reorganization energy, i.e., E'^ = — A. For this particular example, 
this brings the level closer to the chemical potential of the electrodes. As a result, the current 
is enhanced compared with that without vibronic coupling. Although the one-configuration 
result captures this effect qualitatively, it is not quantitative at all, being a factor of three 
too large. 

Figure [2] shows the time-dependent current for larger bias voltages, 0.5V and IV, where 
all other parameters are the same as in Fig. [T] For V = 0.5V, as depicted in Fig. about 
6-8 SPFs per electronic SP group gives the converged result that is in good agreement with 
that obtained with the original model Hamiltonian. The stationary current obtained with 4 
SPFs is 20% larger, which does not satisfy our convergence criterion but is still reasonable. 
The single configuration (1 SPF) result is approximately four times as large as the true 
result. For this bias voltage the artificial correlation effect in the original model becomes 
quite large. It requires 60 SPFs per electronic SP group to obtain the converged result in 
Fig. [2^ for the original model. 

For an even higher bias voltage of = IV, it becomes difficult to converge the steady- 
state current for the original model due to the even increased artificial electronic correlation. 
On the other hand, this also makes the vibrational contribution (relatively) less important. 
Figure Eb shows that within the scattering state representation of the Hamiltonian, the 
steady state current is easily converged, although the transient behavior is different with 
different number of SPFs. In this case, the single configuration result is very accurate, thus 



the steady-state current is very well approximated by the purely electronic model. This is 
due to the rather large voltage. It is known that in the limit of large voltages (within the 
wide-band-approximation), the steady-state current becomes independent of the vibronic 
coupling.*^ For the voltage considered here, the polaron shifted energy level [E^ — A ~ 0.25 
eV) is already well within the bias window given by the chemical potential of the two 
electrodes. 



C. Vibrationally coupled electron transport for Ed = Ef 

Here we consider a model where the energy of the bridge state is located at the Fermi 
energy of the leads. This parameter regime is particularly interesting, because already for 
small bias voltage the transport mechanism corresponds to resonant tunneling and involves 
mixed electron/hole transport. 

In this parameter regime and for a relatively large reorganization energy (e.g., A = 
2000 cm^^), it is found that when transforming to the effective Hamiltonian in the scatter- 
ing state representation, the naive procedure outlined in Sec. IIII A] sometimes gives negative 
steady-state currents, which is unphysical. This numerical problem persists even with a 
relatively large number of scattering channels and a large configuration space. Employing 
the modification described in Sec. IIII Bt on the other hand, gives excellent performance. 
As shown in Figure |3] for a typical range of bath reorganization energies, 1000-2000 cm~^, 
results obtained with 10 SPFs per electronic SP group are in good agreement with those 
obtained with the original model Hamiltonian, which requires 40 SPFs per electronic SP 
group. 

The results depicted in Figs. [3l H] for different values of voltage and electron- vibrational 
coupling strength, all show the same qualitative behavior: a strong increase of the cur- 
rent for short time, which is followed by a decrease to a very small steady-state current. 



As was discussed in detail in in Ref. 



59 



this suppression of the steady state current is a 
manifestation of the phonon blockade effect, which is particularly pronounced for larger 
electronic-vibrational coupling. 

The phonon blockade of the current is usually, qualitatively rationalized by considering 
the energy level of the bridge state. For any finite bias voltage, the bare energy of the bridge 
state {Ed — Ef = 0) is located between the chemical potential of the leads and thus, within 



a purely electronic model, current can flow. The coupling to the vibrations results in a 
polaron shift of the energy of the bridge state. For electronic-vibrational coupling strengths 
A > \V\I2 the polaron-shifted energy of the bridge state is below the chemical potentials of 
both leads and thus the current is blocked. 

The comparison between the converged results and those obtained employing only a 
single SPF for the electronic degrees of freedom (which agree with the results for the purely 
electronic model), demonstrate that this simple picture based on the static energy shift is 
at best qualitatively correct. It is seen from Fig. |3] that the 1 SPF (single configuration) 
approach, i.e., the Hartree-Fock limit in the scattering state representation, significantly 
overestimates the steady state current. The discrepancy becomes larger for stronger electron- 
nuclear couplings. 

It is noted that the results obtained with 1 SPF, which neglect true correlation between 
electronic and vibrational degrees of freedom, are very similar to the results obtained with 
a NEGF method^^i^^ (see, e.g., the data shown in Fig. 9 of Ref. l59l ). In this method, a fac- 
torization of the electronic and vibrational Green's function is employed. As a consequence, 
correlations between electronic and vibrational degrees of freedom are only treated on a very 
approximate level. This results, in this case, in a significant overestimation of the current, 
similar to what would be obtained with a very small (unconverged) number of configurations 
in the ML-MCTDH-SQR method. 

Similar to the examples discussed in Sec. IV Bt the single configuration limit becomes a 
better approximation for a higher bias voltage when the vibronic coupling is fixed. For the 
parameters displayed in Fig. |3l the single configuration result becomes very good for, e.g., 
a bias voltage of 1 V (data not shown), similar to that illustrated in Fig. |2b. Again, in this 
regime, the NEGF approach or the Landauer formula with the shifted bridge state energy 
becomes accurate. 

It should be emphasized, though, that this also depends on the relative strength between 
the vibronic coupling and the bias voltage. Only for voltages \V\ ^ \2{Ed ± A — -E/)|, can 
it be expected that vibronic effects, and thus also vibrationally induced correlation, are 
negligible in the steady state current.— As an example. Figure shows the simulation for 
a higher voltage of 0.5 V but also a larger reorganization energy 4000 cm~^. In this case 
12 SPFs per electronic SP group are required to obtain the result in the scattering state 
representation, where the single configuration result is approximately a factor of three too 



large. Convergence in the original model Hamiltonian also becomes more expensive, which 
requires 64 SPFs per electronic SP group. 

D. Vibrationally coupled electron transport with Ed < Ef 

We finally consider in Fig. [5] an example with a bridge state located at Ed — Ef = —0.5 
eV., i.e. well below the Fermi energy. All other parameters are the same as in Fig. [TJ 
For the voltage considered, this case corresponds to non-resonant transport. Similar to 
the case Ed — Ef = 0, when transforming to the effective Hamiltonian in the scattering 
state representation, the naive procedure outlined in Sec. IIII Al gives a negative steady- 
state current. On the other hand. Figure E] shows the result obtained with employing the 
modification described in Sec. 1111 Bt which gives excellent agreement with the result obtained 
with the original model Hamiltonian. The former requires 8 SPFs to converge, whereas the 
latter require 40 SPFs. The performance with other voltage/vibronic coupling is similar to 
the previous figures. 

VI. CONCLUDING REMARKS 

In this paper, we have shown that the scattering state representation of the underlying 
electronic Hamiltonian represents a promising approach for studying vibrationally coupled 
electron transport through molecular junctions. Within the ML-MCTDH-SQR theory^i^ 
much less configurations are required to achieve convergence than when using the original 
representation of the Hamiltonian. This is because the transformation to the scattering state 
representation effectively diagonalizes the underlying electronic Hamiltonian and, therefore, 
much artificial correlation present in the original model is removed. We have also shown 
that the definition of the scattering states is not unique. This can be used to optimize the 
simulation. Our numerical studies suggest that the form of the Hamiltonian outlined in 
Sec. IIIIBI performs best in all the physical regimes investigated. 

The results of our studies also give physical insight into the importance of correlation 
effects in vibrationally coupled nonequilibrium transport. While for small voltages and/or 
transient effects, electronic-vibrational correlation is often of utmost importance, for higher 
voltages correlation effects are less pronounced in the steady state current. This finding 



is also of relevance for the validity of NEGF calculations for vibrationally coupled electron 
transport, in particular, methods that employ an effective factorization between electronic 
and vibrational degrees of freedom.— 

To conclude this paper, we discuss the computational costs of the simulations in the two 
different representations. As was shown above, much less SPFs are needed when using the 
scattering state representation as compared to the original model. However, the transfor- 
mation used to generate the number operator d~^d and the current operator J^, Eq. f l3.8p . 
scales quadratically versus the number of the electronic states. This is in contrast to the 
original model where the number of electronic operators scales linearly versus the number 
of electronic states. As a result, there are many more operators to evaluate in the scat- 
tering state Hamiltonian, which makes the simulation more expensive. For the examples 
considered in this paper, using 10 SPFs per electronic SP group in the scattering model 
has a similar computational cost to using 40 SPFs in the original model. Therefore, when 
studying vibrationally coupled electron transport at a lower bias voltage, the original model 
Hamiltonian may be more efficient. 

However, things are quite different when the bias voltage becomes higher, where signifi- 
cant artificial correlation exists in the original model, which makes convergence difficult or 
even impossible. The scattering state representation, on the other hand, requires typically a 
smaller configuration space to converge as demonstrated by the numerical examples in this 
paper. This makes the scattering state representation a valuable tool for multiconfigura- 
tional study of current-voltage characteristics where the bias voltage spans a larger range. 
Such work is currently in progress. 

Our current implementation of the scattering state representation has not yet exploited 
its full potential. It is possible that a different energy shift from that of the reorganization 
energy is more efficient for solving the problem in some physical regimes. This may be 
achieved by designing a self-consistent procedure in the convergence tests. Furthermore, a 
different initial condition for the vibrational bath or even a correlated electronic-vibrational 
initial density matrix may be helpful in reducing the transient peaks of the current, and 
thus making the calculation of the steady-state current more stable numerically. Work on 
these issues is in progress. 
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FIG. 1: (a) Time-dependent current I{t) for different representations of tlie model Hamiltonian 
(scattering state vs. original representation), (b) Convergence of I{t) versus the number of SPFs 
for the electronic part, where 6 SPFs are used for each nuclear SP group. The model parameters 
are: = 0.2eV, f3e = leV, and — Ef = 0.5eV for the electronic part; A = 2000cm~^ and 
(jJc = 500cm~^ for the vibrational bath. The source-drain voltage is y = O.IV. 
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FIG. 2: Convergence of I{t) versus the number of SPFs. The parameters are the same as Fig. [T] 
except for the bias voltage: (a) V = 0.5V, (b) V = IV. 
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FIG. 3: Time-dependent current I{t) versus the number of SPFs for the source-drain voltage 
V = O.IV. Six SPFs are used for each nuclear SP group. The model parameters are: = 0.2eV, 
/3e = leV, and E^i — Ej = Q for the electronic part; Wc = 500cm^^ for the vibrational bath and (a) 
A = 2000cm-\ (b) A = lOOOcm-^ 
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FIG. 4: Time-dependent current I{t) versus the number of SPFs for the source-drain vohage 
V = 0.5V. Six SPFs are used for each nuclear SP group. The model parameters are: = 0.2eV, 
/3e = leV, and Ed — Ej = for the electronic part; lOc = 500cm~^ for the vibrational bath and 
A = 4000cm- ^ 
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FIG. 5: Time-dependent current I{t) versus the number of SPFs for the source-drain vohage 
V = O.IV. Six SPFs are used for each nuclear SP group. The model parameters are: Ue = 0.2eV, 
Pe = leV, and E^ — Ej = — 0.5eV for the electronic part; A = 2000cm^^ and ujc = 500cm~^ for the 
vibrational bath. 
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